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Abstract 

In this paper, a real-space analytical expression for the free Green's function (propagator) of 
bilayer graphene is derived based on the effective-mass approximation. Green's function displays 
highly spatial anisotropy with three-fold rotational symmetry. The calculated local density of states 
(LDOS) of a perfect bilayer graphene produces the main features of the observed scanning tunneling 
microscopy (STM) images of graphite at low bias voltage. Some predicted features of the LDOS can 
be verified by STM measurements. In addition, we also calculate the LDOS of bilayer graphene with 
vacancies by using the multiple-scattering theory (scatterings are localized around the vacancy of 
bilayer graphene) . We observe that the interference patterns are determined mainly by the intrinsic 
properties of the propagator and the symmetry of the vacancies. 

PACS numbers: 73.61Wp, 61.72.Ji, 68.37.Ef 
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I. INTRODUCTION 



Since Novoselov et al. fabricated ultra-thin monolayer graphite devices,- the electronic 
properties of a graphite monolayer (graphene) have attracted a great deal of research interest 
due to the Dirac-type spectrum of charge carriers in its gap-less semiconductor material. 
Many interesting properties of single-layer graphene, such as the Landau quantization, the 
defect-induced localization, and spin current states, have been studied experimentally and 
theoretically by several research groups. -^^i^ These non-orthodox properties, including 
massless Dirac fermions around the Dirac points in the first Brillouin zone, result from 
graphene's particular band structure. 

Recent research attention has focused on the electronic structure of bilayer and multilayer 
graphene £.£.iMdLi22.dlM Studies have shown that a bilayer graphene has some unexpected 
properties.— For example, a bilayer graphene shows anomalies in its integer quantum Hall 
effect and in its minimal conductivity on the order of e 2 /h. The common and distinctive 
electronic features of single-layer and bilayer graphene are highlighted in Ref.[15]. Charge 
carriers in a bilayer graphene are mainly quasiparticles with a finite density of states at zero 
energy and they behave similar to conventional non-relativistic electrons. Like the relativistic 
particles or quasiparticles in single-layer graphene, we can describe these quasiparticles by 
using spinor wavefunctions. Although these 'massive chiral fermions' do not exist in the field 
theory, their existence in condensed-matter physics offers a unique opportunity to investigate 
the importance of chirality and to solve the relativistic tunneling problem. 

The unusual physical properties of bilayer graphene are attributed to two key factors, 
(i) The relatively weak inter-layer coupling. Bilayer graphene inherits some properties from 
single-layer graphene material, such as the existence of Dirac points in the first Brillouin 
zone and the degeneracy of electron and hole bands, (ii) The special geometry of bilayer 
graphene with the Bernal stacking (A-B stacking) between adjacent graphene layers. There 
are two kinds of nonequivalent sites (A and B) in each layer as shown in Fig. 1(a) and 1(b). 
Experimental scanning tunneling microscopy (STM) graphite images at low bias voltage 
have verified that only site B is visible and exhibits a triangular structure. In addition, the 
orbital overlap coupling between two adjacent layers is contributed mainly by the carbon 
orbitals at site A.— To date, few analytical studies that attempt to uncover the unique 
electronic properties of the bilayer graphene have been done.— 
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In this paper, we first develop an analytical formula of electronic structure in a bilayer 
graphene with the Bernal stacking based on real-space free Green's function (propagator) 
and the effective-mass approximation. We observe that the physical properties of the bilayer 
graphene are closely related to its propagator and the bilayer graphene behaves similar to a 
massive chiral Fermions system. Since the local density of states (LDOS) can be measured 
by STM, we then compute the LDOS of graphene in various forms based on Green's function 
explicitly. Finally, we highlight the impact of vacancies' patterns (in terms of their geometry 
and symmetry of the computed LDOS) on interferences in STM pictures. 

II. A REAL-SPACE GREEN'S FUNCTION AND ELECTRONIC STRUCTURE OF 
PERFECT BILAYER GRAPHENE 

In what follows, we derive the analytical expression of the free Green's function for bilayer 
graphene in real space based on the effective-mass approximation. Bilayer graphene consists 
of two hexagonal lattice layers coupled by the Bernal stacking as shown in Fig.l (a) and (b). 
In each layer, there are two nonequivalent sites, A and B. Type-A atoms have direct neighbors 
in their adjacent layer, but type-B atoms do not and locate at hollow sites. We assume that 
there is one valence electron per carbon atom in a bilayer graphene. Provided that the 
difference between wavefunctions orthogonal to the graphene plane can be neglected, the 
bilayer graphene can be modeled as an effective two-dimensional system. Since the physical 
properties of graphene are determined by tt bands near the Dirac point, only the contribution 
from the n band is considered in this paper. We further limit our analysis by considering 
only the nearest interactions of the graphene's p z orbitals. The tight-binding Hamiltonian 
of the bilayer graphene is 



where (v\i)i, I = (1,2) is a wavefunction at site i at layer /. 62 P is the on-site energy, V\ 
is the nearest hopping parameter within the layer and V2 is the nearest hopping parameter 




+^(|*)llO'| + |*>220'|+/*.C.) 




(1) 
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between two layers. Our presented calculations are conducted with e 2p = OeV, V\ = S.OeV 
and V 2 = OAeV. 

The Bloch orbits for two nonequivalent sites, A and B, are written as \kA)i = 
-jjf Yl e ik ' r ^ \j A )i , \k B )i = -j= ^e 4k ' rj s \j B )i. The summation covers all sites A and B on 

3A 3B 

layer I. ta and r# denote the real-space coordinates of sites A and B, respectively. Here, N 
is the number of unit cells in the crystal. The Hamiltonian can be rewritten as 

I e 2p Vifi* V 2 \ 

VijJ, e 2p 

V 2 e 2p Vxn 

\ Vm* e 2p ) 



H 



(2) 



where /i = e lkya + e %<K ~ 



■{ \Hiak x a ky ^ g 

+ e H 2 —> and a = 1.42A is the c-c bond length. By 
defining the retarded Green's function as, G r Q et = l/(e — H + irj), where 77 is an infinitesimal 
quantity, we obtain the reciprocal space Green's function for bilayer graphene 



G r et \ 



( t(t 2 -V?nn*) V^if - V?(i(i*) 
V^{t 2 - V?w*) t(t 2 - V 2 ^* - V 2 ) 
t 2 V 2 



tv 1 v 2 ^* 

V 2 V 2 fi* 2 



t 2 V 2 W1V2II \ 

WMn v 2 v 2 ^ 2 
t(t 2 - vfnn*) Vm(t 2 - vfnn*) 
Vxfi*(t 2 - V?w*) t(t 2 - v?w* - v 2 2 ) J 

(3) 



where A = (t 2 - Vfn/i* - tV 2 )(t 2 - \ffifj,* + tV 2 ) and t = e - e 2p + irj. From A = 0, the 
dispersion relation is expressed as e = ±-y ± J (7A;) 2 + (-^) 2 with 7 = 3aVi/2, which is 
consistent with the previous theoretical result.— 

The real-space Green's function of the bilayer graphene can be derived by taking the 
Fourier transform of (jQ e *(k, e). For simplicity, these calculations are carried out around 
six corners in the first Brillouin zone within the low-energy region based on the effective- 
mass approximation,. Six Dirac points can be divided into two equivalent sets of points, 
K X ,K 3 ,K 5 and K 2 ,K 4 ,K 6 , which are shown in Fig. 1(c). They form two 360° integrals 

around K 1 and K 4 . By using the mathematical relation 

00 

e ik.(r M -0 = Mklr ^ _ r /j) + 2 j>j B (A ; |p„ - r'J)cos(n^ K ), 

n=l 

where (/i, v) = (A, B), <p r r '„ is the angle between k and r M — t' v , and J n is the type-I n-order 
Bessel function. We can simplify the real-space Green's function of the top layer (I = 1) to 



G r et (v A ,v' A ,e) = cos[K 1 -(r A -r' A )]-F 1 (\r A -v' A \,e) 
G r et (r B ,r' B ,e) = cosfK 1 • (r B - r' B )}-F 2 (\r B - r' B \, e) 
G r et (r A , r' B , e) = smjK 1 • (r A - r' B ) + a rA yJ 

■F 3 (\r A -r' B \,e) 
G r et {r B , r' A , e) = sin\K} ■ (r B - r' A ) - a rB yJ 

■F 3 (\r B -r' A \,e), (4) 

Here K 1 = (4^/3^,0), cosifr^y = cos(6 — ot r r ' v ), 9 is the angle between the k and x 
axis. cx r tr > v is the angle between r M — r' u and x axis. The definitions of Fi,F 2 , and F 3 are 



F 1 (\r,-T , v \,e) = - f C dkkJ {k\v,-v' u \) 
71 Jo 



t 

+ 



H 2 - ^k) 2 - tV 2 t 2 - ( 1 k) 2 + tV 2 i 
F 2 {\v ti -r'Xe) = -[ dkkJ (k\v^-v' u \) 



t-V 2 t + V 2 



t 2 -{ ri k) 2 -W 2 t 2 -(-fk) 2 + tV 2 

k-c 



F 3 (|r M -^|,e) = -/ Wjx^lr^-r'J) 

JO 



7 7 



h 2 - ink) 2 - tv 2 t 2 - ( 7 £;) 2 + tv 2 r 



(5) 



where S = 3^/3a 2 /2 is the area of the unit cell in real space and k c is the cutoff wave 
vector. The corresponding cutoff wave length 2n/k c is comparable to the lattice constant 
a. With the same approach in our previous work to solve electronic structure in the single- 
layer graphene, k c is set to be 1.71/a.— From Eq.(4), it is clear that the real-space Green's 
function of the bilayer graphene is constructed by multiplying two terms. The first term is 
spatially anisotropic, which can be represented by sine and cosine functions determined from 
two nonequivalent sets of the Dirac points as shown in Fig. 1(c). The second term including 
Fx, F 2 or F 3 is spatially isotropic in real space and depends only on the distance between 
two sites. 
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Many physical properties of bilayer graphene can be deduced from this explicit expression. 
Interestingly, when k c — ► oo, the functions, F\, F 2 and F 3 , have simple analytical forms 

F^-r'^e) = — 2 [tK,C-\v,-v' u \) 

7T7 7 

+tK C-^\r,-v'M 
F 2 (|r,-r'J, £ ) = ^[(t-V^-lr.-r'J) 

7T7 Z 7 

+ (t + y 2 )K (^|r M -r'J)], 
7T7 7 

+i* 2 ifi(— |r M - r^D], (6) 
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where t x = v^ 2 - ^ and t 2 = v^ 2 + tV 2 . K and K x are the zero-th order and the 
first order terms of the type-I modified Bessel function, respectively. Compared with our 
previous study— the expressions of the Green's function derived here are the same as those 
of the single-layer graphene once the nearest hopping parameter between two layers (V2) is 
set to zero. That is to say, it is straightforward to compare theoretical results of bilayer 
graphene with previous theoretical and experimental data of the monolayer graphene by 
setting V 2 = 0M 

Fig. 2 shows the LDOS of bilayer graphene with e=— 0.1 eV. Here, the LDOS of site r M 
is calculated by setting Pq{y^e) = — -Im(jo e *(r M , r M , e). We can clearly observe the contrast 
between the LDOS at site A and that at site B. Site B is highlighted and forms a triangular 
lattice with three-fold symmetry. According to the Tersoff-Hamann model^ 1 which has been 
successfully used to explain experimental results,— 1 ^ the STM image can be simulated using 
the LDOS of the sample surface. It is reasonable for us to compare the LDOS of bilayer 
graphene with previous experimental STM observations of graphite surfaces at low bias 
voltage (i.e. 0.3 V). Our calculated LDOS can clearly produce the main features as shown 
in several observed STM images reported in Ref. [16-18]. 

One direct approach to find the difference between bilayer and monolayer graphene is 
to study the LDOS at sites A and B for both monolayer and bilayer graphenes. Fig. 3(a) 
and 3(b) present the LDOS of one carbon site (A or B) in the monolayer graphene and 
two sites (A and B) in the bilayer graphene, respectively. For a monolayer graphene, sites 
A and B are equivalent and thus only one curve of LDOS is shown in Fig. 2(a). However, 
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these two sites are nonequivalent for a bilayer graphene. Fig. 2(b), therefore, clearly shows 
that the LDOS of bilayer graphene at site B is greater than that at site A, when |e| is less 
than approximately 0.4 eV. This LDOS difference leads to the brighter spots at site-B as 
shown in Fig. 2. When |e| > 0.4 eV, the difference between the LDOS at sites A and B in the 
bilayer graphene gradually diminishes. This feature is easy to verify through STM dl/dV 
mapping at a relatively large bias voltage. For example, sites A and B can be both observed 
in STM dl/dV mapping with the applied bias voltage of 0.6 V. This phenomenon can also 
be elaborated by our analytic expression of Green's function. 

The difference between the LDOS of site A and that of site B can be defined as 



After some simple math operations, shown in Appendix Eq.(A-5) and Eq.(A-6), Ap (e) can 
be simplified to 



From this equation, the LDOS of the bilayer graphene can be divided into two regions 
based on the interlayer hopping parameter One region is |e| < V2, where the LDOS of 
site B is larger than that of site A by a constant, Tjpp-. The other region is \e\ > V2, where the 
LDOS of site A equals that of site B. Our calculations agree well with our numerical results 
when V2 = 0.4 eV. The calculations show the dispersion relation is approximately quadratic 
dependent in the low-energy region and the results match the previous theoretical results.— 
We find that the low-energy LDOS of both sites A and B are linearly proportional to the 
energy (e) (refer to the Appendix). Note that the LDOS of one unit cell at the very low- 
energy level, however, is a constant. The LDOS can be expressed as p(e) = [po( r yi, £ ) + 
Po{tBi Vbi £)] — which is the same as those in the conventional two-dimensional system. 
This remarkable feature of the bilayer graphene offers a direct method to measure the inter- 
layer hopping parameter V2, which is also the threshold needed to distinguish LDOS contrast 
in STM images. The intra-layer hopping parameter V\ can be deduced based on the LDOS 
difference measured according to STM images ( Ap = ^% with 7 = 3aVi/2). 



Apo(£) 



-Im[G r et (r B , t b , e) - G r et (r A , r A , e)] 

7T 

-Im[F 2 {0,e)-F 1 {0,e)). 



(7) 
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III. LOCAL DENSITY OF STATES OF BILAYER GRAPHENE WITH LATTICE 
VACANCY 



To validate the accuracy of real-space Green's function derived in this paper, we use 
the promoted formula to solve the electronic structure of a bilayer graphene with a single 
vacancy. Based on the tight-binding scheme, a vacancy can be simulated by introducing a 
large on-site energy at the vacancy site on a bilayer graphene.— Assuming that the single 
vacancy locates at site B with on-site energy U, the transport matrix or the T-matrix— ^ 
can be written as T = U(l — UG™ 1 )' 1 . If U is large, the T-matrix becomes 



T(0,0,e) ~ -G^' 1 (0,0, e) 

= -F^fre). (9) 

Here, the position of the vacancy is set as the coordinate origin. That is to say, Rb = in 
the x — y plane. Using the Dyson equation, we have 

G ret (Tfj,, r„, e) = GQ et (r fJj ,r u ,e) 

+G r et (r„0,e)T(0,0,e)G r et (0,r„,E), 

(10) 

where G ?re *(r /Lt , r u , e) is the Green's function of bilayer graphene with a single vacancy in real 
space. The LDOS on site r M can be determined by 

p(r„,e) = —ImG^^r^e), 

7T 

which can be simply rewritten as 

p(r A ,e) = p (r A ,e) + sm 2 (K 1 • r A + a rA )H x {x A ,e), 

p(r B ,e) = p (r B ,e) + cos 2 ^ 1 ■ r B )H 2 (r B ,e). (11) 



In the above equation, 



po{r A ,e) = ~Jm[Fi(0,e)], 

7T 

po{r B ,e) = Im[F 2 (0,e)}, 

71 

H ^ - \ Im ^f v (12) 



Fig. 4(a) shows the calculated LDOS with a single B-site vacancy on the top sheet of the 
bilayer graphene. An interesting feature is clearly observed; bright spots localized near the 
vacancy around site A have nice three-fold rotational symmetry. This phenomenon can be 
explained based on Eq.lfTll. Since the function H 2 in the Eq.fTHi) has a very small value 
(close to zero), its magnitude changes little as the distance from the vacancy increases and 
thus the magnitude of p(rs,£) remains very small. Although the cosine term in Eq. (TTll 
reaches its maximum when tb = ^^n, n = 0, ±1, ±2 • ■ - , the sine function in p(r^,e) 
reaches its maximum along directions with angles, a rA = 30°, 90°, 150°, 210°, 270° and 330°. 
Compared to H 2 , the value of the function Hi has a relatively larger value, i.e. 0.4 eV~ x at 
l.o A away from the vacancy for e = -0.1 eV. The function Hi, however, decays rapidly with 
increasing distance from the vacancy. Note that sites A along the directions a TA = 90°, 210° 
and 330° are located closer to the vacancy than those along the directions a TA = 30°, 150° 
and 270°. The calculated site- A LDOS value is between that of the nearest and that of the 
next nearest to the B-site vacancy, that is 0.56 eV -1 and 0.36 eV -1 , respectively. Therefore, 
the sites along the a TA = 90°, 210° and 330° directions are brighter than those along the 
a rA = 30°, 150° and 270° directions. These bright spots show the localized character of 
the region surrounding the vacancy. By comparing to the dl/dV image of a vacancy on 
the single-layer graphene^, the enhanced localization of LDOS in the bilayer graphene is 
caused by the additional inter-layer channel propagated from the vacancy point. Green's 
function can oscillate in a significantly longer distance starting from the vacancy in the 
monolayer graphene than in the bilayer graphene. The asymptotic behavior of Hi can help 
us understand this phenomenon. When the distance is longer than 4.0 A, Hi is less than 0.1 
eV -1 and thus these bright spots show the localized character of the region within nearest 
lattice surrounding the vacancy. 

The LDOS with vacancy near site- A is also simulated as shown in Fig. 4(b). We observe 
the similar features that the vacancy around site-B has bright spots with three-fold symme- 
try. The calculated site-B LDOS value is between that of the nearest and that of the next 
nearest to the site-A vacancy, which is between 0.35 eV -1 and 0.23 eV -1 . 

The extension of our calculation to consider several vacancies on bilayer graphene is 
straightforward. The scattering T-matrix in Eq.(10) includes all contributions from va- 
cancies. Similar to a single-vacancy case, a large on-site energy U for all vacancy sites is 
used in our simulations. Fig. 4(c) shows the LDOS near an AB-type vacancy (a couple of 
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the neighboring sites A and B) on a bilayer graphene surface. The LDOS is still localized 
around the vacancies. However, the clear three-fold rotational symmetry disappears in this 
case. The bright spots localized above the vacancy correspond to the B sites. Note that 
the LDOS at the bright spots is about 0.03 eV -1 , which is less than those in Fig. 4(a) and 
(b) by an order of magnitude. This phenomenon results in the destructive interference due 
to the multiple-scattering caused by two vacancies. The asymmetric pattern in Fig. 4(c) 
reflects the residual contribution from the destructive interference since two sites A and B 
are nonequivalent in the bilayer graphene. Our result suggests that the interference pattern 
is determined mainly by vacancies, provided that the spatial anisotropy of the Green's func- 
tion is fixed. The symmetry will be lost in LDOS whenever vacancies break the three-fold 
rotational symmetry in bilayer graphene. 



IV. CONCLUSION 

In this paper, an analytical form of the real-space Green's function (propagator) of bi- 
layer graphene is constructed based on the effective-mass approximation. Green's function 
demonstrates an elegant spatial anisotropy with a three-fold symmetry for defect-free bilayer 
graphene. The LDOS of the bilayer graphene determines the main features of experimental 
STM images on graphite surfaces with low-bias voltage. The predicted features according to 
our simulated results can be verified by STM measurements. For example, two nonequivalent 
atomic sites can be observed in STM dl/dV images with different bias voltages. The infor- 
mation of interlayer and intralayer hopping strength can be deduced based on the contrast 
of STM images. Moreover, We also calculate the LDOS of bilayer graphene with vacancies 
by using the multiple-scattering theory. The interference patterns are determined mainly 
by the properties of Green's function and the symmetry of the vacancies. Once the vacan- 
cies break the intrinsic symmetry of the graphene, the three-fold rotational symmetry of the 
LDOS vanishes. Our model provides exact results near the Dirac points. We have discovered 
some interesting STM patterns of the second layer, but these results will be published else- 
where. In this paper, the bilayer graphene is described by using the simple non-interactive 
tight-binding scheme. We are currently investigating how Coulomb interaction impacts the 
electronic structure of the bilayer graphene. 
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APPENDIX 



From Eq.(3), Green's functions at the top layer in reciprocal space are expressed as 



2H 2 - Vf/i/i* -tV 2 t 2 -Vffifi* +tV 2 



2 



2 l t 2 - V?nn* -tV 2 t 2 - V 2 fifx* + tV 2 



G oab(^ £ ) - o(t^ t/2....* 777" + 



Tt 2 - Vffifi* -tV 2 t 2 - Vffifi* + W 2 1 



G 0Ba(^ £ ) - -[-g T/2 ....* 777 + 



Tt 2 - Vf/i/i* -tV 2 t 2 - Vf/i/i* + tV 2 1 



(A-l) 



By taking the Fourier transform of G^ik, e) in the first Brillouin zone (1BZ), we can obtain 
the exact expression of Green's function in real space for bilayer graphene Go 6 *(r M , r' v , e) = 
f 1BZ dkGn^(k, e)e l ^ r ^~ r ' v \ where \x and v are for site-A or site-B atoms, respectively. Based 
on the effective-mass approximation, we can sum k points near the six corners (labeled by 
i — 1 ~ 6) in the first Brillouin zone. The real-space bilayer graphene Green's function can 
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be written as 

cn , 4 , r ; 4 , £) = _g_ ± J dKd ^ _ m + fi _ _ _ ]e ^>^ 

= E / ^1 P - . tV , + P - ^ + tVa le^--.' 

(A-2) 

In Eq.(A-2), the integral around if 1 , K 3 , if 5 and if 2 , if 4 , if 6 can be summed together 
separately to form two 360° integrals around if 1 and if 4 , that is 



f27T 

x / ^' t2 _ (7 ; )2 _ ft/2 + f2 _ (7 ; )2+tv . i e "-'-^ 



^0 



x 



dd f C dkk\ l ~ ^ + t + V2 l e * k -(r B -r' s ) 

J I k[ t 2 -( 1 k) 2 -tv 2 + t 2 -( 1 k) 2 + tvJ e 



G r Q et (v A ,v' B ,e) = jA-A^^ 



(2tt) 2 

x d9l dkk 2 -y(-isin9 - cos6)[- — — L - — + - — — L - U^^'b) 

1 1 H t 2 - (7A;) 2 - tV 2 t 2 - (7A;) 2 + tv 2 



jo 

if 



+e iK 4 (r A -r' s ) 
»2tt r k, 

x 1 



q 

r<rett I \ _ a ( JK 1 (T B -T' i ) 

^0 { r B,rA, £ ) - J^y \ e 



x del dkk^iisinO - cos6)[- — — — — + -, — r-rr-, ]^B-r' A ) 

Jo Jo ,v t 2 -{ 1 k) 2 -W 2 t 2 -(jk) 2 + tV 2 l 



+e iK*(r B -r' A ) / dQ 

Jo 

x £ dkkHismO + cose )[ t2 _ {i l )2 _ tV2 + t2 _ {l l y + tv / ^} 

(A-3) 
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By using relation 

oo 



n=l 



and integrating the angle part of Eq.(A-3), we can easily get Eqs.(4) and (5). If we let 
k c — > oo, Eq.(5) is simplified to Eq.(6). 

Next, we briefly derive Eq.(8). At the same site A or B, the corresponding Green's 
functions are: 

r-iretf \ ^ /^jnr t — V 2 t + V 2 , 



(A-4) 



The LDOS at sites A and B are 

p (r A , r A ,e) = --Im[G r et (r A , r A , e)} 

71 



S j/mjf e d[^( 7 fc) 2 + (y) 2 ] 



2(tt 7 ) 



t t 

{ 1 , + 

e + i V -f-J (7A;) 2 + (f ) 2 e + ztj - f + ,/( 7 A;) 2 + (f ) 2 



2 V V /'V 1 V 2 / ' 2 1 V V r" ) 1 V 2 



£ + IT] + f - \J (7A;) 2 + (-y) 2 £ + i77 + f + ^/( 7 A:) 2 + (f ) 2 



} 



S|e| 

7T7 2 

2tt7 2 
1 



\e\ > V 2 

. (A-5) 

e < Vo 



p (r B ,r B ,e) = Im[G r e (r B , r B , e)] 

7T 



2(tt 7 ) 



t-v 2 t-v 2 

+ 



'e + iv-^ - ^(7^ + (f ) 2 ^ + ^-f + v/(7A:) 2 + (f ) 2 
t + V 2 t + F 2 



£ + «77 + f - y/( 7 A;) 2 + (f ) 2 £ + MJ + f + ^( 7 fc)2 + (^)2 

lei > ^ 2 



} 



S|e| 

7T7 2 " 

g(|e|+X2) 

27T7' 



3— |£ < V2 



(A-6) 
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Then, the difference between LDOS at A and B sites is 



Apo(e) = Po(r B ,r B ,e) - p (r A ,r A ,e) 
f Id > V 2 




(A-7) 
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Figure 1: (color online) The crystal structure of a bilayer graphene. The unit cell consists of 
two layers with two nonequivalent sites: A (yellow) and B (gray), (a) Top view with surface 
unit vectors al and a2, (b) Side view, (c) The first Brillouin zone of a bilayer graphene, where 
K 1 ~K 6 are the Dirac points. These Dirac points can be further divided into two nonequivalent 
sets, K 1 ~ (K\K 3 ,K 5 ) and K 4 ~ (K 2 ,K 4 ,K 6 ). 
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Figure 3: (a) LDOS of monolayer graphene on one site (sites A and B are equivalent), (b) LDOS 
of bilayer graphene on two sites (sites A and B are nonequivalent). 
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Figure 4: (color online) LDOS of vacancies in bilayer graphene. (a) single B-site vacancy; (b) A-site 
vacancy; (c) a pair of AB vacancy. Here,e=-0.1 eV 
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